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1.0  SUMMARY 


Recent  experiments,  in  which  voids  are  precision  cut  using  a  second  laser,  have  been  used  to 
study  void  coalescence.  In  that  study,  both  the  Rice  and  Tracey  model  and  the  Thomason  model 
were  shown  to  have  reasonable  agreement  with  experiment  results  for  the  copper,  in  most  cases. 
However,  no  simulations  were  performed  in  which  voids  were  explicitly  represented.  This  article 
uses  finite  element  simulations  to  describe  void  coalescence  as  observed  in  copper  bar  uniaxial 
tension  experiments.  Several  problems  occur  when  meshing  the  three-dimensional  geometry 
containing  176  voids,  namely  element  distortion  and  damage,  minimum  time  step,  and 
appropriate  material  model  parameters.  Further,  a  temporal  and  spatial  convergence  study  was 
used  to  estimate  errors,  thus,  this  study  helps  to  provide  guidelines  for  modeling  of  materials 
with  voids.  Finally,  we  use  a  Gurson  model  with  Johnson-Cook  strength  to  simulate  the  void 
growth.  Simulations,  using  the  codes  ABAQUS  explicit,  EPIC,  and  Presto,  agree  well  with 
experiments. 

Keywords:  Coalescence;  Fracture;  Simulations;  Damage  Modeling 

2.0  INTRODUCTION 

Damage,  in  solid  mechanics,  is  the  state  of  a  material  which  has  lost  structural  integrity  or  other 
mechanical  properties  of  interest.  Material  damage  can  involve  void  coalescence.  This  may  be 
especially  true  for  ductile  fracture  mechanisms  (see  for  example  ref.  [1]).  Ductile  type  fracture 
behavior  may  occur  when  shear  bands  form,  especially  at  elevated  strain  rates  as  observed  in  the 
Kalthoff  experiment  [2],  which  is  an  experiment  that  is  challenging  to  simulate  due  to  the 
difficulty  of  correctly  modeling  damage. 

Damage  can  involve  void  nucleation,  growth,  and  coalescence,  which  are  physical  phenomena 
that  require  high  resolution  computational  approaches  for  simulation.  Experiments  used  to  study 
these  physical  phenomena  include  gas-gun,  expanding  ring  fragmentation,  penetration,  and 
perforation  or  plugging.  The  range  of  length  scales  needed  to  span  void  nucleation  and  crack 
formation  are  on  the  order  of  104  (micrometers  to  centimeters).  Cracks  that  have  formed  may 
grow,  bifurcate,  merge  with  other  cracks,  or  arrest. 

To  visualize  void  growth  and  coalescence  in  a  controlled  manner,  Week  et  al.  [3]  experimented 
with  model  materials  that  contain  laser-drilled  holes.  Void  growth  and  coalescence  were 
observed  via  X-ray  tomography  with  impressive  resolution.  This  behavior  is  important  when 
attempting  to  understand  damage  evolution.  The  tensile  test  that  was  used  is  sufficient  to 
examine  stress  triaxiality  and  local  plastic  behavior,  which  both  help  to  isolate  the  relevant 
physics  for  a  relatively  simple  experiment.  Therefore,  this  experiment  was  selected  as  an 
appropriate  benchmark  problem. 

We  are  motivated  to  simulate  void  growth  and  coalescence  to  validate  codes  1)  to  better  validate 
continuum  void  growth  models,  2)  to  know  what  computational  strategy  will  best  simulate  the 
physics,  and  3)  to  better  understand  the  limitations  of  our  codes.  The  nearly  pure  copper  material 
is  used  in  this  article  because  we  have  Johnson-Cook  model  parameters.  Moreover,  pure  copper 
is  a  ductile  metal  with  many  applications  in  which  damage  occurs. 
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3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 


3.1  Void  Coalescence 

Void  coalescence  has  been  modeled  at  the  macroscale  continuum  level  by  several  authors  [4], 

[5],  [6],  [7],  [8],  [9],  Tonks  et  al.  incorporated  void  linking  via  an  intervoid  ligament  as  a  basic 
mechanism  that  was  used  to  model  the  local  instability  [4],  Kim  et  al.  [5]  studied  the  effects  of 
stress  triaxiality  on  void  coalescence  using  the  Gurson-Tvergaard  constitutive  relation.  Another 
Gurston-type  approach  was  proposed  by  Lassance  et  al.  [6]  and  used  to  study  penny-shaped 
voids  resulting  from  particle  fracture.  Void  nucleation,  growth  and  coalescence  were  modeled 
by  Hammi  and  Horstemeyer  [7]  using  an  anisotropic  damage  progression  formulation,  in  which 
coalescence  was  modeled  as  a  second  rank  tensor  that  is  governed  by  the  plastic  strain  rate  tensor 
and  the  stress  state,  where  the  coalescence  threshold  is  related  to  the  void  length  scale.  Maire  et 
al.  [8]  modeled  their  uniaxial  tension  experiments  using  the  classical  Rice  and  Tracey  approach 

[10]  and  later  Landron  et  al.  [9]  used  the  Rice  and  Tracey  approach  that  was  corrected  by  Huang 

[11] ,  In  the  later  article,  cavity  shape  change  was  incorporated  and  used  to  improve  model  fits  to 
experiment  results. 

Finite  element  simulations  of  void  coalescence  were  performed  by  some  of  the  aforementioned 
authors.  Kim  et  al.  [5]  presented  a  finite  element  approach  that  included  a  representative  unit 
cell  containing  a  spherical  void  at  its  center.  Lassance  et  al.  [6]  studied  a  three  dimensional 
periodic  packing  of  voids  by  a  cylinder  with  appropriate  boundary  conditions.  In  that  study,  void 
coalescence  was  possible  before  necking  in  a  purely  uniaxial  tension  condition,  provided  that  the 
initial  void  volume  fraction  was  sufficiently  high.  Several  other  authors  have  simulated  a  single 
pore  geometry  using  finite  elements.  In  one  of  the  more  recent  studies,  Siad  et  al.  [12]  studied 
void  coalescence,  but  with  a  sole  focus  on  finite  element  simulations  of  a  single  cylindrical  void 
geometry.  Ha  and  Kim  [13]  investigated  void  growth  and  coalescence  in  an  anisotropic 
crystalline  material,  namely  f.c.c.  single  crystals,  using  the  finite  element  method.  In  their  study, 
the  growth  and  coalescence  behavior  were  simulated  using  a  single  pore. 

Previous  void  coalescence  experiments  have  been  conducted.  See  for  example,  the  non¬ 
destructive  X-ray  tomography  experiments  in  refs.  [3],  [8],  [9].  Recent  experiments  by  Week  [3] 
isolate  void  growth  and  coalescence  in  a  controlled  manner  via  laser-drilled  holes.  In  contrast, 
natural  void  evolution  and  coalescence  was  observed  by  Landron  et  al.  [9],  where  damage 
evolution  in  dual-phase  steel  was  studied  [8],  [9],  In  the  latter  study,  the  correction  by  Huang 
[1 1]  to  the  Rice  and  Tracey  model  [10]  was  validated  for  different  triaxiality  states. 

3.2  Code  Verification 

Verification  is  the  process  of  determining  that  a  model  implementation  accurately  represents  the 
developer’s  conceptual  description  of  the  model  and  the  solution  to  the  model  [14],  Four 
predominant  sources  of  error  are  1)  insufficient  spatial  discretization  convergence,  2)  insufficient 
temporal  discretization  convergence,  3)  lack  of  iterative  convergence,  and  4)  computer 
programming  [15].  Programming  errors  and  iterative  convergence  are  not  addressed  in  this 
paper.  The  solution  over  the  entire  computational  domain,  including  boundaries,  must  be  verified 
for  the  geometry  and  loading  conditions  of  interest. 
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None  of  the  previous  studies  mentioned  here  included  an  explicit  mesh  for  multiple  voids. 
Moreover,  void  coalescence  was  not  simulated  for  an  experimental  geometry  that  included  voids. 
Instead,  the  aforementioned  studies  focused  on  meshing  a  single  void,  and  in  one  case,  periodic 
boundary  conditions  were  used  [6],  These  studies  only  used  one  mesh  and  the  largest  allowable 
time  step,  thus,  neither  spatial  nor  temporal  refinements  were  used.  Therefore,  this  article  is  the 
first  to  attempt  to  quantify  spatial  and  temporal  convergence  rates,  asymptotic  regime,  and 
corresponding  errors.  These  activities  are  part  of  uncertainty  quantification,  which  addresses  the 
three  fundamental  components  of  computer  simulations  for  physical  systems,  namely,  model 
qualification,  model  verification,  and  model  validation  [14].  Although  previous  penetration 
studies  are  strong  in  model  qualification,  model  validation  is  impossible  without  complete 
verification,  i.  e.  numerical  error  quantification. 

The  grid  convergence  index  (GCI)  is  a  relatively  simple  method  to  estimate  errors  due  to  mesh 
refinements.  However,  GCI  methods  cannot  provide  statistical  confidence,  which  can  be 
obtained  using  response  surface  methods  (RSM)  [16].  Another  limitation  for  GCI  methods  is  the 
use  of  an  empirical  safety  factor  Fs  to  provide  a  confidence  interval.  Some  of  the  assumptions 
made  in  the  GCI  and  RSM  methods  are  relaxed  in  a  nonlinear  ansatz  error  model  [17]. 

Confidence  estimation  for  a  given  safety  factor  is  based  on  the  number  of  grid  points  used  Ng 
and  no  consensus  has  been  reached  on  the  value.  Approximately  95%  certainty  (5%  uncertainty, 
that  is  roughly  a  2cr  error  band  if  the  distribution  is  Gaussian)  is  assumed  for  Ng  >  3,  and 
Fs  =  1.25  [18].  Although  these  values  are  for  steady  state  fluid  flow  and  heat  transfer,  they  are 
commonly  used  in  a  wide  range  of  hydrocode  verification  studies.  A  consensus  has  been 
reached  that  multiple  methods  should  be  explored  in  a  systematic  verification  study  and  that  GCI 
methods  ( Ng  >  3)  often  produce  useful  information.  For  example,  relative  GCI  values  for 
different  grid  sets  were  used  in  Ref.  [19]  to  determine  the  best  grid  set  based  on  finding  the 
lowest  GCI. 

3.3  Verification  Approach 

Verification  involves  quantifying  the  error  associated  with  solving  the  governing  equations 
regardless  of  the  values  of  the  model  parameters.  The  most  important  part  of  verification  is  to 
perform  spatial  and  temporal  refinement  studies. 

The  first  step  is  to  identify  the  metrics  or  quantities  of  interest  for  observing  convergence, 
namely  1)  displacement  ( Ac )  and  2)  load  or  force  ( Lc ),  at  the  end  of  the  bar  at  which  point  voids 

coalesce.  Experiments  described  above  provide  accelerometer  data  and  penetration  depth. 
Therefore,  penetration  depth  is  considered  to  be  the  most  important  verification  metric  in  this 
study. 

3.3.1  Models  of  Convergence  Error. 

The  exact  solution  is  denoted  as  F*  and  discrete  solutions  are  denoted  in  simplified  notation, 

Fx  t  =  F (Ax,  At).  We  define  the  general  error  metric  in  Eq.  (1)  using  the  norm  of  the  difference 
between  the  continuous  and  numerical  solutions: 
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For  the  metrics  used  in  this  paper,  the  error  is  simply  the  absolute  value  of  the  difference,  that  is, 
ex,t  =  \F*  ~  Fx  t | . 

Discretized  solutions  that  neglect  higher  order  terms  according  to  the  Richardson  extrapolation 
estimation  (REE)  technique  are  given  by, 

FXlt  =  F*  +  a*hf  (2) 

where  a*  is  a  fitting  constant,  p*  is  the  convergence  rate  or  order  of  convergence,  and  is  either 
the  spatial  or  temporal  step  size,  1  /NceUs.  Equation  (2)  is  rarely  if  ever  observed  in  practice. 
Alternatively,  the  approximate  solution  F  approximates  F*  by  extrapolation  such  as  the  REE 
technique  or  by  response  surface  methods.  Therefore,  we  can  approximate  Eq.  (2)  by, 

Fx,t  =F  +  ahf  (3) 

where  the  REE  fitting  parameters  are  {1,  a,  p).  The  corresponding  space-only  or  time-only  error 
estimation  is  given  by, 


ex,t  =  2o  +  dhf  +  h.  o.  t.  (4) 

Higher  order  terms  ( h .  o.  t.)  are  neglected  in  this  three  parameter  fit  with  parameters  {!0,  a,  p). 

The  asymptotic  region  is  described  in  this  study  to  be  the  set  of  spatial  and  temporal  mesh  sizes 
that  provide  a  consistent  and  accurate  order  of  convergence  and  GCI.  This  region  is  generally 
unknown  before  performing  simulations.  An  estimated  order  of  convergence  may  be  calculated 
using  three  grid  points  Ng  =  3,  constant  refinement  ratio  rh,  constant  a,  and  constant 
convergence  rate  p,  as  in  Roache  [20],  as 


Pr  =  log 


Or  ~  F3) 

(Fi  -  F2) 


/log  [rh] 


(5) 


where  F1,  F2,  and  F3  are  the  fine,  medium,  and  coarse  grid  solutions,  respectively.  If  the  grid 
convergence  is  monotonic  with  constant  convergence  rate,  then  pr  =  p*  and  the  extrapolated 
solution  is  given  by, 


F  =  F*  =  Ft  + 


(Ei  +  F2) 

~  1) 


(6) 


However,  F  =  F*  is  generally  not  true.  Further,  p  and  a  are  rarely  constant.  A  safety  factor  Fs  is 
used  to  provide  the  GCI,  which  is  expressed  as 


GCI  =  Fq 


E.  —  F7 


Fi 


/(>f  - 1) 


(7) 
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A  conservative  safety  factor  Fs  =  3  is  used  with  Ng  =  3  is  assumed  to  estimate  the  uncertainty 
68%  of  the  time  (or  lcr)  because  we  assume  a  Gaussian  error  distribution. 

Three  main  assumptions  are  made  in  the  formulation  above.  The  error  according  to  Eq.  (4) 
assumes  that,  1)  numerical  solution  convergence  is  monotonic,  2)  space-time  coupling  effects  are 
neglected,  and  3)  higher  order  terms  are  negligible. 

Coefficients  in  Eq.  (4)  are  solved  analytically  as  in  Eq.  (5).  Throughout  this  paper,  we  use  the 
following  error  definition: 


&x,t  —  F  Fxt 


(8) 


3.3.2  Grid  Convergence  Study. 

Several  uniaxial  tension  simulations  are  used  to  investigate  spatial  and  temporal  refinements. 
The  following  two  steps  are  used  for  the  nominal  case  refinement  study: 

Extrapolate  solution  F  from  Eq.  (6)  assuming  p  =  pr. 

Apply  the  three  parameter  space-only  fit  using  Eq.  (4)  and  optimization  procedure  described 
above. 


The  spatial  refinement  set  is  dx  =  [W/10,  W/20/,W/40],  where  W  is  the  width  of  the  tension 
specimen.  The  temporal  refinement  set  is  CFL  =  [1.0, .5, .25],  where  the  Courant-Friedrichs- 
Lewy  (CFL)  is  a  fraction  of  the  maximum  stable  time  step,  or  in  this  article,  the  fixed  time  step 
in  which  mass-scaling  was  used.  Although  the  CFL  does  not  correspond  directly  to  the  time 
step,  CFL  values  are  appropriate  for  temporal  refinements  in  tensile  loading  simulations  because 
the  minimum  time  step  is  determined  from  conditions  near  the  edges  of  the  pores.  Therefore, 
pore  growth  and  coalescence  are  sensitive  to  temporal  refinements  using  CFL  values. 

The  stable  time  step  using  ABAQUS  was  on  the  order  of  le-1 1  seconds,  which  is  prohibitively 
small  even  on  a  high  performance  computing  cluster.  Therefore,  we  used  mass-scaling. 


3.4  Macroscopic  Modeling 

In  addition  to  the  simulation,  in  which  the  mesh  is  refined  such  that  geometry  of  the  voids  is 
explicitly  included,  we  compare  a  macroscale  damage  model  approach.  For  the  EPIC  code,  we 
use  the  following  Gurson-type  yield  criterion  (see  for  example  Tvergaard  and  Needleman  [21]). 


o>  = 


f 

+ 2 qx  cosh 

V 


2cr 


M  j 


-l-<h2=0 


(9) 


where  the  macroscopic  Mises  stress  is  <Je  =  (3  SyS9  /2j  ,  in  terms  of  the  stress  deviator 

Sy  =  criJ  -  Gij o\  / 3 ,  and  o\  is  the  macroscopic  mean  stress.  For  qx  =  1 ,  Eq.  (9)  is  that  derived 

by  Gurson  [1], 
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A  Mie-Gruneisen  equation  of  state  (EOS)  is  used.  The  pressure  for  the  virgin  material 
[22]  is  given  by, 


P  =  P(p,e)  =  PH 


f 


1 


A 


1 - Yp  +  Ype 

v  2  J 


(10) 


where  p  =  p/(l-</>)  is  the  density  of  the  solid  material,  <j)  is  the  ratio  of  void  volume  over  the 
total  volume,  e  =  e  is  the  internal  energy  of  the  solid  material,  p  -  p  I  p()  - 1  is  the  current 
density  p  and  the  initial  density  p(),  T  is  the  Gruneisen  constant  for  the  solid  material.  The 
Hugoniot  pressure  consistent  with  the  standard  quadratic  EOS  form  is  given  in  Eq.  (1 1)  in  terms 
of  three  parameters  that  are  fit  to  experiments,  namely  kx ,  k2 ,  and  A\ . 


PH  —  k\P  +  k2p^  +  k2p^  (1 1) 

For  EPIC,  Presto,  and  Abaqus  simulations,  we  used  the  Johnson-Cook  strength  model  [22], 
which  is  given  by, 


<r  =  [A  +  Benp\L  +  C\n.e  *][l-T*m]+czP  (12) 

where  s  is  the  equivalent  plastic  strain,  s*  =  s/s0  is  the  dimensionless  total  strain  rate  where 

e  is  the  equivalent  total  strain  rate,  s0  =  1  .Os-1 ,  T*  =  (T  -  Twom  )j(jmelt  -  Troom  )  is  the 

homologous  temperature,  P  is  the  hydrostatic  pressure,  and  A,  B,  n,  C,  m,  and  a  are  material 
constants  [22], 

3.5  Finite  Element  Simulations 

The  tensile  test  that  was  used  in  Ref.  [3]  is  used  here  as  an  appropriate  benchmark  problem.  The 
full  tensile  specimen  is  dimensioned  in  Figure  1  (a).  The  gage  length  is  the  portion  of  the 
specimen  used  in  the  simulations  and  is  the  2.0  mm  length  shown  in  Figure  1  (b).  The  specimen 
with  a  gauge  length  rectangular  cross  section  was  constructed  from  three  plates.  The  10pm 
middle  plate  had  176  holes  with  10pm  diameters  drilled  using  a  laser  in  a  rectangular  array 
oriented  at  45°  to  the  tensile  axis.  The  center-to-center  hole  spacing  was  20um,  and  further 
details  of  the  array  are  shown  in  Figure  2. 

The  simulated  boundary  conditions  are  built-in  ends  that  load  the  bar  in  uniaxial  tension  at  a 
strain  rate  of  0.1  s'1. 
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Figure  1:  Uniaxial  Tension  Specimen  (a)  Actual  Specimen;  (b)  Gage  Length 


Figure  2:  Laser  Drilled  Rectangular  Array  of  176  holes  (only  partially  shown) 
Oriented  at  45°  to  the  Tensile  Axis 
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The  hydrocodes  Abaqus  [23],  Presto  [24],  and  EPIC  [25]  are  used  in  the  current  work  for  their 
leading  capabilities.  These  codes  were  used  with  Lagrangian  implementation,  appropriate 
constitutive  models,  and  equation  of  state.  Additionally,  all  three  codes  are  capable  of  using  a 
macroscale  continuum  model  similar  to  Gurson’s  model. 

Abaqus  and  Presto  both  have  approximately  second  order  accuracy  in  space  and  time. 
Lagrangian  methods  in  both  codes  are  used  with  tetrahedron  elements  in  this  study.  Further 
details  can  be  found  in  Ref.  [26]  In  an  effort  to  minimize  computational  time,  the  specimen 
geometry  was  initially  meshed  with  hexagonal  elements;  however,  it  was  found  to  be  not 
possible  to  further  refine  the  hexagonal  mesh  in  the  porous  regions.  The  mesh  was  therefore  re¬ 
generated  with  tetrahedral  elements  in  the  porous  region  of  the  specimen  only,  leaving  the  un- 
porous  top  and  bottom  regions  meshed  with  hexagonal  elements.  This  scheme  was  also  found 
problematic  due  to  the  use  of  a  tie  constraint  by  the  hydrocodes  at  the  transitional  surfaces  in 
order  to  combine  the  different  element  regions  in  the  model.  These  tie  constraints  could 
introduce  artificial  constraints  to  the  model,  which  could  lead  to  a  restraining  behavior  in  the 
specimen  from  necking.  It  was  decided  to  maintain  the  most  realistic  integrity  of  the  model  and 
accept  the  longer  computational  time  by  generating  a  full  tetrahedron  mesh  of  the  entire 
specimen. 

For  all  simulations,  we  used  the  Johnson-Cook  (HJC)  model  [27]  with  Mie-Griineisen  equation 
of  state  (EOS)  to  model  the  pressure,  density,  and  internal  energy  relationship.  Element  failure 
was  assumed  to  occur  when  i)  elements  inverted,  ii)  nodal  Jacobian  ratio  became  negative,  or  iii) 
the  element  equivalent  plastic  strain  exceeded  150%. 

In  Abaqus,  the  equivalent  plastic  strain  rate  is  defined  as, 


n 


(13) 


for  a  >  <7° 


where  a  is  yield  stress  at  nonzero  strain  rate,  a°(spl,  0,f{)  is  static  yield  stress,  9  is  temperature, 
fi  are  other  field  variables,  and  D(0,  /[)  and  n(0,fi)  are  material  parameters. 

Material  failure  or  element  deletion  is  mesh  dependent.  Although  material  failure  is  mesh 
dependent,  solution  convergence,  if  it  exists,  is  not  necessarily  the  same  as  for  the  conservation 
equations  and  momentum  balance.  Errors  can  be  reduced  by  refining  the  mesh  in  locations  in 
which  material  is  expected  to  fail.  Therefore,  refinements  were  made  to  all  elements  in  the 
simulated  specimen,  specifically  those  in  the  porous  region  where  elements  tend  to  fail  first.  The 
specimen  geometry  was  partitioned  into  three  regions  as  shown  in  Figure  3  below.  The  mesh 
refinement  was  established  by  decreasing  the  global  mesh  seed  of  the  entire  model.  Additional 
refinement  was  performed  on  Region  2  of  the  specimen,  where  the  refinement  metric  used  was 
the  number  of  element  line  segments  per  circular  void.  Table  1  below  lists  the  refinement  details 
for  all  three  models  generated. 
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Figure  3:  Uniaxial  Tension  Specimen  (Model  3) 


Table  1:  Mesh  Refinement  Summary 


Region 

Model  1  (fine) 

Model  2  (middle) 

Model  3  (course) 

Nodes  Elements 

Nodes 

Elements 

Nodes 

Elements 

1 

28,620 

152,359 

13,076 

65,813 

6,078 

28,036 

2 

107,957 

621,838 

47,061 

267,045 

21,980 

122,093 

3 

28,669 

152,574 

13,130 

66,173 

6,038 

27,775 

Totals 

165,246 

926,771 

73,267 

399,031 

34,096 

167,904 

Void 

Refinement 

12  Element  Line  Segments 

8  Element  Line  Segments 

6  Element  Line  Segments 
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Figure  4:  Uniaxial  Tension  Specimen  of  the  Coarsest  Mesh  (Model  3).  (a)  Full  Model  of 
Specimen;  (b)  Close  up  of  the  Cross-Section  of  the  Porous  Region 

The  macroscale  continuum  void  region  is  estimated  based  on  the  6.5%  void  volume  fraction 
specified  by  Week  [3].  This  region  has  length  and  width  approximated  based  on  the  size  of 
rectangular  void  region  shown  in  Figure  5  above.  The  width  and  height  of  a  rectangle  tangential 
to  the  outside  rows  of  voids  is  160  microns  and  307  microns,  respectively.  A  buffer  region  was 
placed  around  this  tangential  rectangle  to  create  an  approximate  continuum  void  region  with 
width  and  height  of  170  microns  and  311  microns,  respectively.  The  actual  voids  had  a  through 
thickness  of  10  microns;  however,  to  maintain  the  total  volume  of  these  voids  for  a  6.5%  void 
volume  fraction,  the  thickness  of  this  continuum  void  region  was  set  to  46.3  microns. 


Units  in 
microns 


— 170- 

-160 

Vs— 


Vs- 


307  311 


Figure  5:  Dimensions  of  Continuum  Void  Region 
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4.0  RESULTS  AND  DISCUSSIONS 

The  following  equations  from  Ref.  [3]  were  used  to  generate  results  in  all  simulations: 


a  =  s(e  +  1) 

(14) 

s  =  ln(e  +  1) 

(15) 

L 

a  =  — 

(16) 

A 

s-mfe) 

(17) 

where  Eqs.  (14)  and  (15)  are  used  to  calculate  true  stress-strain  curves  prior  to  necking,  and  Eqs. 
(16)  and  (17)  are  used  to  calculate  these  curves  after  necking.  In  these  equations,  s  =  L/A0  is 
the  engineering  stress,  e  =  (l  —  /0)//0  is  the  engineering  strain,  L  is  the  load,  A0  is  the  initial 
cross-sectional  area,  and  h,  and  l  are  the  initial  and  current  gage  lengths,  respectively.  For  the 
true  stress  vs.  true  strain  calculations,  the  area  in  the  middle  of  the  specimen  where  the  necking 
occurred  was  calculated  by  tracking  the  node  displacements  at  the  comers  during  each 
simulation.  This  tracked  area  was  used  in  Eqs.  (16)  and  (17)  to  compute  true  stress  and  true 
strain  values. 

Results  for  the  two  metrics  are  given  for  various  mesh  refinements  in  Table  2.  From  Ref.  [3],  the 
true  strain  at  necking,  sn  =  0.26,  crn  =  250  MPa,  is  found  from  the  intersection  between  the  true 
stress  and  the  work-hardening  rate  curves  as  a  function  of  true  strain,  as  shown  in  Figure  6. 

We  obtain  for  each  metric,  a)  order  of  convergence,  and  b)  extrapolated  solution  with 
corresponding  uncertainty,  and  c)  computation  time  for  each  simulation.  Order  of  convergence 
and  GCI  are  compared  to  results  obtained  from  experiments.  Convergence  results  and  the 
solution  uncertainty  calculated  for  CFL  values  are  summarized  in  Table  3. 

Stress  vs.  strain  curves  from  FEA  simulations  are  shown  in  Figure  7,  Figure  8,  and  Figure  9  for 
comparison  with  the  experimental  stress  vs.  strain  curves  shown  in  Figure  6.  Note  that  we 
compare  only  the  Copper  experimental  data,  not  Glidcop.  ABAQUS  simulations  for  Models  1, 

2,  and  3  are  shown  in  Figure  7  (a),  (b),  and  (c),  respectively.  Presto  simulation  Models  2  and  3 
are  shown  in  Figure  8  (a)  and  (b),  respectively.  Model  1  simulations  for  Presto  were  never 
completed  due  to  a  very  large  estimated  run  time  as  shown  in  Table  4,  which  was  predicted  after 
15%  of  completion  of  the  simulation.  EPIC  simulations  for  Models  1,  2,  and  3  are  shown  in 
Figure  9  (a),  (b),  and  (c),  respectively. 


Table  2:  Mesh  Dependency  Summary  (OFHC  Copper) 


dx 

Abaqus 

EPIC 

Presto 

sn  an  (MPa) 

<?n  (MPa) 

<7„(MPa) 

Model  1 

0.2262 

319.74 

0.2104 

310.81 

- 

- 

Model  2 

0.2271 

320.64 

0.2083 

310.20 

0.2143 

314.14 

Model  3 

0.2282 

321.39 

0.2030 

308.58 

0.2030 

314.20 
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Table  3:  Convergence  Study,  L  =  low  (rh  =  2.187),  H  =  high  (rh  =  2.329) 


dx 

Abaqus 

EPIC 

<7„  (MPa) 

<7„(MPa) 

Pr  (L/H) 

1.020/0.945 

0.775/0.717 

1.610/1.490 

1.657/1.533 

Fext 

0.2255 

318.7 

0.2112 

311.0 

GCI 

0.0098 

0.0019 

0.0119 

0.0022 

Table  4:  Simulation  Run  Times 


Abaqus 

EPIC 

Presto 

run  time 
(hrs) 

run  time 
(hrs) 

run  time 
(hrs) 

Model  1 

200.8 

401.8 

952.4 

Model  2 

53.3 

204.9 

328.1 

Model  3 

18.6 

109.9 

253.5 

True  strain 


Figure  6:  True  stress  vs.  true  strain  Week’s  experiment. 
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Stress  (MPa)  Stress  (MPa) 


(a) 


(b) 


(c) 


Figure  7:  Stress  vs.  Strain  for  Abaqus  Simulations. 


(a) 


(b) 


Figure  8:  Stress  vs.  Strain  for  Presto  Simulations. 
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Figure  9:  Stress  vs.  Strain  for  EPIC  Simulations. 


Void  coalescence  is  shown  in  Figure  10,  where  void  screenshots  were  taken  at  0%  true  strain  and 
at  50%  true  strain  from  Abaqus  Model  1  simulation  and  are  compared  with  the  tomography 
pictures  from  the  experiment  in  Ref.  [3],  The  stretch  of  the  voids  from  Abaqus  simulation  very 
closely  matches  the  experimental  pictures. 
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o  o  o 
o  o  o  c 

20  nm 


(a) 


20  iJirn 

(c) 


•  •  I  < 

4  •  I 

20pm 

(b) 


20fim 

(d) 


Figure  10:  Abaqus  Void  Visualization  vs.  Experiment  at  0%  True  Strain  in  (a)  and  (b),  and 

at  50%  True  Strain  in  (c)  and  (d) 


5.0  CONCLUSIONS 

OFHC  Copper  material  properties  from  the  EPIC  material  library  were  used  for  all  simulations. 
The  true  stress  at  necking  for  all  simulations  was  found  to  be  around  60-70  MPa  higher  than  the 
experiment  from  Ref.  [3]  as  Figure  6  shows.  This  seems  to  indicate  that  the  simulations  are  not 
capturing  the  elastic  response  of  the  material  accurately  and  the  material  model  used  in  the 
simulations  therefore  does  not  match  the  material  properties  of  the  specimen  tested  in  Ref.  [3], 
In  order  to  test  the  codes  capability  of  predicting  void  coalescence  and  failure  more  accurately, 
the  material  model  should  be  first  calibrated  with  an  experimental  test  with  no  voids  present  in 
the  specimen.  Once  the  material  model  is  accurate,  then  the  computational  simulations  for  the 
specimen  with  voids  will  compare  more  realistically  with  experimental  results. 

It  was  seen  in  this  work  that  the  best  mesh  to  use  this  study  would  be  Model  2.  This  is  based  on 
the  consideration  that  the  Abaqus  simulation  was  0.71%  sn  error  and  0.60%  cn  from  fully 
converged  solutions  based  on  the  mesh  convergence  results,  and  that,  for  the  EPIC  simulations, 
the  errors  were  1.39%  for  en  and  0.26%  for  an.  Use  of  Model  2  is  advised  over  Model  1  since 
the  Abaqus  simulation  run  time  for  Model  2  is  73%  faster  than  for  Model  1  and  the  EPIC 
simulation  the  run  time  for  Model  2  is  49%  faster  than  for  Model  1 . 

The  effects  of  element  deletion  were  not  closely  studied  in  these  simulations;  however,  it  should 
be  mentioned  here  that  the  deletion  criteria  for  the  Abaqus  simulations  was  purely  specified  by 
the  Johnson-Cook  material  failure  model,  while  in  the  Presto  and  EPIC  simulations,  an 
equivalent  plastic  strain  failure  of  300%  was  specified  as  an  additional  deletion  criteria.  As  far 
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as  studying  the  number  of  deleted  elements  with  changing  mesh  refinement,  there  is  not  data 
available  from  these  simulations. 

Note  that  all  three  codes  produced  such  very  similar  results  that  we  cannot  base  a  comparison  of 
prediction  capability  on  these  results  alone.  If  run  times  are  considered  as  a  metric  to  judge  the 
codes’  performance,  then  Abaqus  achieves  the  fastest  run  times  out  of  the  three. 
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